package br.edu.ufcg.msnlab2.grupo07.monteCarlo;

import java.util.ArrayList;
import java.util.HashMap;
import java.util.List;
import java.util.Map;
import java.util.Map.Entry;

import mathExpr.UnknownDefinitionException;
import br.edu.ufcg.msnlab2.MSNLabException;
import br.edu.ufcg.msnlab2.misc.Function;


/**
 * This class is responsible for implementing a version of the Monte Carlo
 * Method for calculating an integral value using the concept of recursive stratified
 * sampling. This algorithm is referred as the MISER Monte Carlo.
 * @author David Candeia
 * @author Felipe Leal
 */
public class MiserMonteCarloSolverImpl extends MonteCarloSolver implements MiserMonteCarloSolver, PlainMonteCarloSolver {

	private final double PFAC = 0.1;
	private final long MNPT = 15;
	private final long MNBS = 60;
	private final double TINY = 1E-30;
	private final double BIG = 1E30;
	
	private long iran = 0;
	
	//Results to be returned
	private double average = 0d;
	private double var = 0d;
	
	/**
	 * This method will request the recursive execution of the MISER algorithm in order to approximate the
	 * integral value
	 * @param function The function expression
	 * @param limits The limits of the integral for each variable
	 * @param numberOfPoints The number of points to be sampled
	 * @return The integral estimative and its error
	 * @throws MSNLabException Exception thrown if any invalid condition occurs while evaluating the function
	 * in the points being considered
	 */
	public MonteCarloResult solve(Function function,	Map<String, Tuple<Double, Double>> limits, long npts, double dith) throws MSNLabException {
		if(dith < 0 || dith > 0.5){
			throw new MSNLabException("Dithering must be between 0 and 0.5");
		}
		
		iran = 0;
		
		MonteCarloResult result = solveMiserRecursively(function, limits, npts, dith);
		
		double volume = calculateVolume(limits);
		result.setIntegrationResult(result.getIntegrationResult() * volume);
		result.setError(result.getError() * volume);
		return result;
	}
	
	/**
	 * This method execute the MISER algorithm in order to approximate the integral value
	 * @param function The function expression
	 * @param limits The limits of the integral for each variable
	 * @param numberOfPoints The number of points to be sampled
	 * @param dith Parameter that causes dimension to not be divided exactly in the middle. For a middle division it should be zero.
	 * @return The integral estimative and its error
	 * @throws MSNLabException thrown if any invalid condition occurs while evaluating the function
	 * in the points being considered
	 */
	public MonteCarloResult solveMiserRecursively(Function function, Map<String, Tuple<Double, Double>> limits, long numberOfPoints, double dith) throws MSNLabException {
		long npre = 0;
		long nptl = 0;
		long nptr = 0;
		double fracl = 0d;
		double s = 0d;
		double sigl = 0d;
		double siglb = 0d;
		double sigr = 0d;
		double sigrb = 0d;
		double sum = 0d;
		double sumb = 0d;
		double rgl = 0d;
		double rgr = 0d;
		double rgm = 0d;
		
		if(numberOfPoints <= 0){
			numberOfPoints = MNBS-1;
		}
		
		int ndim = limits.size();
		
		//Too few points, use plain Monte Carlo
		if(numberOfPoints < MNBS){
			Map<String, List<Double>> sampledPoints = samplePoints(limits, numberOfPoints);
			calculateFunctionValueAndVariance(sampledPoints, function, numberOfPoints);
			
			//Returning result
			MonteCarloResult result = new MonteCarloResult(var, average);
			return result;
		}else{//Inspect for bisection
			List<Double> rmid = new ArrayList<Double>();
			npre = Math.max((long)(numberOfPoints*PFAC), MNPT);//Points to be used in dimension divisions variance estimation
			List<Double> fmaxl = new ArrayList<Double>();
			List<Double> fmaxr = new ArrayList<Double>();
			List<Double> fminl = new ArrayList<Double>();
			List<Double> fminr = new ArrayList<Double>();
			
			//Initializing upper and lower bounds for dimensions
			int index = 0;
			for(Entry<String, Tuple<Double, Double>> entry : limits.entrySet()){
				iran = (iran*2661+36979) % 175000;
				if(iran-87500 < 0){
					s = -dith;
				}else{
					s = dith;
				}
				
				double lower = entry.getValue().getK();
				double upper = entry.getValue().getA();
				double mid = (0.5+s)*lower + (0.5-s)*upper;

				rmid.add(index, mid);
				
				fminl.add(index, BIG);
				fminr.add(index, BIG);
				
				fmaxl.add(index, -BIG);
				fmaxr.add(index, -BIG);
				index++;
			}
			
			//Calculating function values for each sampled point
//			System.out.println("Points to sample "+npre);
			Map<String, List<Double>> sampledPoints = samplePoints(limits, npre);
//			System.out.println("Sampled Points: "+sampledPoints.toString());
			
			for(int i = 0; i < npre; i++){
				double result = function.evaluate(sampledPoints, i);
				
				//Evaluating the lowest and highest function value obtained for each dimension 
				int listIndex = 0;
				for(Entry<String, List<Double>> entry : sampledPoints.entrySet()){
					if(entry.getValue().get(i) <= rmid.get(listIndex)){
						fminl.add(listIndex, Math.min(result,fminl.get(listIndex)));
						fmaxl.add(listIndex, Math.max(result,fmaxl.get(listIndex)));
					}else{
						fminr.add(listIndex, Math.min(result,fminr.get(listIndex)));
						fmaxr.add(listIndex, Math.max(result,fmaxr.get(listIndex)));
					}
					listIndex++;
				}
			}
			
			//Evaluating which dimension has the higher variance
			sumb = BIG;
			int jb = 0;
			siglb = 1d;
			sigrb = 1d;
			
			for(int k = 0; k < ndim; k++){
				Double maxl = fmaxl.get(k);
				Double minl = fminl.get(k);
				Double maxr = fmaxr.get(k);
				Double minr = fminr.get(k);
				
				if(maxl > minl && maxr > minr){
					sigl = Math.max(TINY, Math.pow(maxl-minl, 2d/3d));//Variance approximation in the left
					sigr = Math.max(TINY, Math.pow(maxr-minr, 2d/3d));//Variance approximation in the right
					
					sum = sigl + sigr;
//					System.out.println("Dimension "+k+" variance "+sum);
					if(sum <= sumb){
						sumb = sum;
						jb = k;
						siglb = sigl;
						sigrb = sigr;
					}
				}
			}
			
//			System.out.println("Dimension chosen "+jb);
			
			//Clearing some lists
			fminr.clear();
			fminl.clear();
			fmaxr.clear();
			fmaxl.clear();
			
			//Verifying if any dimension gave a variance value
			if(jb == 0){
				jb = (int)((ndim*iran)/175000);
			}
			
			//Finding limits for dimension to be bisected
			rgm = rmid.get(jb);
			index = 0;
			for(Entry<String, Tuple<Double, Double>> entry : limits.entrySet()){
				if(index == jb){
					rgl = entry.getValue().getK();
					rgr = entry.getValue().getA();
					break;
				}else{
					index++;
					continue;
				}
			}
			
//			System.out.println("Dimesion limits: ("+rgl+", "+rgm+", "+rgr+")");
			
			fracl = Math.abs((rgm-rgl)/(rgr-rgl));
			nptl = (long)( MNPT+(numberOfPoints-npre-2*MNPT)*fracl*siglb/(fracl*siglb+(1-fracl)*sigrb) );
			nptr = numberOfPoints - npre - nptl;
			
			//Preparing limits for left recursion
			Map<String, Tuple<Double, Double>> leftRecursionLimits = new HashMap<String, Tuple<Double,Double>>();//Map that will contain limits for recursion
			leftRecursionLimits.putAll(limits);
			index = 0;
			for(Entry<String, Tuple<Double, Double>> entry : leftRecursionLimits.entrySet()){
				if(index == jb){
					Tuple<Double, Double> newTuple = new Tuple<Double, Double>(entry.getValue().getK(), rgm);
					entry.setValue(newTuple);
					break;
				}
				index++;
			}
			
//			System.out.println("Left recursion with limits "+leftRecursionLimits.toString()+" and "+nptl+" points ");
			MonteCarloResult leftResult = solveMiserRecursively(function, leftRecursionLimits, nptl, dith);
			
			//Preparing limits for right recursion
			Map<String, Tuple<Double, Double>> rightRecursionLimits = new HashMap<String, Tuple<Double,Double>>();//Map that will contain limits for recursion
			rightRecursionLimits.putAll(limits);
			index = 0;
			for(Entry<String, Tuple<Double, Double>> entry : rightRecursionLimits.entrySet()){
				if(index == jb){
					Tuple<Double, Double> newTuple = new Tuple<Double, Double>(rgm, entry.getValue().getA());
					entry.setValue(newTuple);
					break;
				}
				index++;
			}
			
//			System.out.println("Right recursion with limits "+rightRecursionLimits.toString()+" and "+nptr+" points ");
			MonteCarloResult rightResult = solveMiserRecursively(function, rightRecursionLimits, nptr, dith);

			//Composing results
			double integrationResult = fracl*leftResult.getIntegrationResult()+(1-fracl)*rightResult.getIntegrationResult();
			double error = Math.pow(fracl, 2)*leftResult.getError()+Math.pow((1-fracl),2)*rightResult.getError();
//			System.out.println("Combinend estimative: "+integrationResult+" cmbined error estimative "+error);
			
			return new MonteCarloResult(error, integrationResult);
		}
	}
	
	/**
	 * This method execute the MISER algorithm in order to approximate the integral value
	 * @param function The function expression
	 * @param limits The limits of the integral for each variable
	 * @param numberOfPoints The number of points to be sampled
	 * @return The integral estimative and its error
	 * @throws MSNLabException thrown if any invalid condition occurs while evaluating the function
	 * in the points being considered
	 */
	public MonteCarloResult solve(Function function,	Map<String, Tuple<Double, Double>> limits, long numberOfPoints) throws MSNLabException {
		return this.solve(function, limits, numberOfPoints, 0);
	}
	
	/**
	 * This method is responsible for calculating the function value and variance used in the Miser implementation.
	 * @param sampledPoints The set of sampled points.
	 * @param function The function being evaluated
	 * @param numberOfPoints The number of sampled points
	 * @throws MSNLabException thrown if any invalid condition occurs while evaluating the function
	 * in the points being considered
	 */
	private void calculateFunctionValueAndVariance(Map<String, List<Double>> sampledPoints, Function function, long numberOfPoints) throws MSNLabException {
		
		if(sampledPoints.size() == 0){
			throw new RuntimeException("Invalid set of sampled points!");
		}
		
		for(String variable : sampledPoints.keySet()){
			if(sampledPoints.get(variable).size() == 0){
				throw new RuntimeException("Invalid set of sampled points!");
			}
		}
		
		if(numberOfPoints <= 0){
			throw new MSNLabException("Invalid number of points!");
		}
		
		if(function == null){
			throw new RuntimeException("Invalid function!");
		}
		
		try{
			int index = 0;
			long maxIndex = numberOfPoints;

			double summ = 0d;
			double summ2 = 0d;
			var = 0d;
			average = 0d;
			
			
			for(int i = index; i < maxIndex; i++){
				
				//Calculating function value for a point
				double result = function.evaluate(sampledPoints, i);
				summ += result;
				summ2 += Math.pow(result,2);
			}
			
			average = summ / numberOfPoints;
			var = Math.max(TINY, (summ2 - summ*summ/numberOfPoints)/Math.pow(numberOfPoints, 2));  
		
		}catch(UnknownDefinitionException e){
			throw new RuntimeException(e.getMessage());
		}
	}
}
